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Abstract 

An improved method is presented for the numerical evaluation of multi-loop inte- 
grals in dimensional regularization. The technique is based on Mellin-Barnes repre- 
sentations, which have been used earlier to develop algorithms for the extraction of 
ultraviolet and infrared divergencies. The coefficients of these singularities and the 
non-singular part can be integrated numerically. However, the numerical integration 
often does not converge for diagrams with massive propagators and physical branch 
cuts. In this work, several steps are proposed which substantially improve the behavior 
of the numerical integrals. The efficacy of the method is demonstrated by calculating 
several two-loop examples, some of which have not been known before. 



1 Introduction 



Perturbation theory is one of the most important tools for the calculation of particle physics 
observables. In many cases, radiative corrections are important to match the experimental 
precision. General and efficient techniques for the computation of one-loop corrections have 
been developed over many years, resulting in automated frameworks for the calculation of 
processes with up to six external legs [TJ. Most methods for one- loop corrections reduce 
the loop integrals to a well-defined set of master integrals by Passarino-Veltman and Mel- 
rose reduction [2], integration-by-parts and Lorentz-invariance relations [3], or on-shell cut 
techniques pEj. See Ref. [5] and references therein for a more detailed overview. 

In some cases, however, one-loop corrections are not sufficient and two- or even three- 
loop corrections becomes relevant. Typical examples are electroweak precision observables 
[HIE] and production of gauge bosons [8], top quarks [HI EI] or Higgs bosons [TT] at the 
Large Hadron Collider (LHC). Many highly complex calculations have been performed for 
this purpose by several authors (see references in previous sentence), involving specialized 
techniques tailored to each particular process. For two-loop amplitudes with a small number 
of different mass and momentum scales one can follow the traditional approach of reduction 
to a small number of master integrals, which are then calculated analytically. However, this 
strategy does not work for applications with a large number of different scales, due to the 
complexity of the expressions and the fact that in general the master integrals cannot be 
determined analytically beyond the one-loop level. 

Fully numerical techniques for the evaluation of two- and higher-loop integrals face two 
main challenges: extraction of ultraviolet (UV) and infrared and collinear (IR) singularities, 
as well as stability and efficiency of the numerical integrations. Methods based on dispersion 
relations [H[l2] and integrations over Feynman parameters [13] result in low- dimensional 
numerical integrals with good convergence properties, but they do not offer a general al- 
gorithm for the isolation of divergencies. Two powerful approaches for a straightforward 
extraction of 1/(4 — D) poles in dimensional regularization are sector decomposition [H] 
and Mellin-Barnes representations [15141?] . However, these methods lead to complicated 
multi-dimensional numerical integrals, which have convergence problems for diagrams with 
physical branch cuts. For sector decomposition, the instability related to physical branch 
cuts can be circumvented by the integration contour deformation proposed in Ref. [TS] . 

This article instead focuses on the numerical integration of Mellin-Barnes (MB) represen- 
tations. Compared to sector decomposition, MB representations of Feynman integrals have 
the advantage that they typically lead to substantially shorter expressions in the integrand. 
However, MB representation will only become a viable general tool for numerical two-loop 
integrations if the behavior of the numerical integrals can be substantially improved. In this 
paper, several techniques for this purpose are presented, which can also be used in combina- 
tion to optimize the convergence properties of the integrals. In particular, complex variable 
transformations and relations that reduce the number of integration dimensions prove to be 
very useful. 

Section [2] starts out by reviewing the employment of MB representations for loop integrals 
and explaining the origin of the numerical instabilities in the MB integrals. In section [3] the 
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new methods to deal with these instabilities are described in detail. Numerical results for 
several examples of two-loop integrals are presented in section HI Where applicable, they are 
compared to existing results in the literature, which have been obtained by other methods. In 
addition, results are shown for some integrals appearing in the calculation of two-loop QCD 
corrections to gg — > ti which have not been calculated before. Finally, the main findings are 
summarized in section [5j 



2 Mellin- Barnes representations 



The following discussion of the construction of MB representations and the isolation of UV 
and IR singularities is largely based on Ref . [16] . Let us start with a one-loop n-point integral 
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where M = [Iq + • • • + \i n . The integrand denominator can be transformed into a MB 
representation by 
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where the following two conditions must be met: 



1. The integration contours Ci for Zi are straight lines parallel to the imaginary axis chosen 
such that all arguments of the gamma functions have positive real parts. 

2. The Ai are complex numbers with | arg(Aj) — arg(A,-)| < ir for any pair i, j. 

Satisfying the second condition potentially requires some adjustments of the Feynman parametriza- 
tion. For example, the scalar diagram in Fig. [1] is represented by the Feynman integral 
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Figure 1: Scalar box diagram corresponding to eq. (pi}. Thick lines indicate propagators 
with mass m, while thin lines are massless, and external legs are constrained to be on-shell. 

where the terms with m 2 violate condition [2] above. More generally, problems with the 
application of the MB representation arise if one mass or momentum parameter appears 
both with positive and negative signs in the denominator of the Feynman integral. This 
complication can be avoided by using the relation Y2i x i = 1- ^ n the case of eq. 01]), this 
relation allows to rewrite it as 

/f^= T(4 — y ) [ dxodx\dx 2 dxz — (5) 

Jo [— SX1X3 — txoX2 + m?XQX2 + m 2 x\ — ie] ~ 

so that all terms with m 2 now have positive signs. The term —SX1X3 does not lead to any 
difficulties since one can imagine to temporarily shift s slightly off the real axis into the 
complex plane so that \arg(— sxiXs)\ < tt is satisfied. Thus the Mellin-Barnes formula can 
be safely applied to eq. (jSJ). 

After the introduction of the MB representation, the integration over the Feynman pa- 
rameters can be carried out straightforwardly, using 

C dx --- dx n 6(1 - x - . . . - x n ) x^- 1 ■ ■ ■ c- 1 = r M---rK) (6) 

Jo r(ct + ... + a n ) 

assuming that the exponents satisfy 

Re(oi) > 0. (7) 

A MB representation for a two-loop integral can be obtained by first deriving a MB 
representation for one sub-loop, as described above, and inserting the result into the second 
loop, which then can be transformed into a MB representation with the same procedure. In 
general, the MB representation for the first sub-loop will contain invariants depending on 
the second loop momentum, which can be interpreted as new propagators for the second 
integration. 

As an example, consider the "sunset" diagram in Fig. [2J The MB representation for the 
upper sub-loop is given by 

1^ = j^y J dz,dz 2 (ml)-^-**(mir(-q 2 2 r r(- Zl )r(-z 2 )T(l + z 1 + z 2 ) 

Hi - e - z 2 )T(e + z x + zo) 
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Figure 2: "Sunset" two-loop diagram. As before, all propagators are assumed to be scalars. 



where e = (4 — D)/2 and q 2 is the internal momentum of the second (lower) loop. For the 
second integration one then obtains 

'" m= pS? / dzidz '< I [ [?2r " [(92 +pf ~ ml 
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so that the {q 2 ) zi gives rise to a new propagator. After constructing the MB representation 
for the second loop, the final result is 

$ m = f ^i^^K)-^-^^) 22 ^) 1 -^ 1 - 2 ^-/) 23 r(-z 2 )r(-* 3 ) 

X T(l + + z 2 )T{z 3 - Zx) — —— . (10 

F(2 - e + z 3 ) 

The MB integral can then be be performed numerically. 

As mentioned above, the Mellin-Barnes formula and the Feynman parameter integration 
require the arguments of all gamma functions to have a positive real part. In general, this 
requires e to differ from zero by a finite amount. When taking the limit e — > 0, the poles 
of some gamma functions can cross some of the integration contours, so that one needs to 
correct for this by picking up the residue of the integrand for crossed each pole. As a result, 
for e — > one ends up with the original integral plus a sum of terms with lower integration 
dimension stemming from these residues. An algorithm for this procedure is described in 
Refs. [HUE]. 

Since for appropriate choices of the contours each integral in this sum is regular, the 1/e 
singularities must be contained in the coefficients of these integrals. Therefore this method 
allows to extract all UV and IR divergencies in a systematic way. 

3 Improving numerical convergence 

The contour integrations in the MB representation can be parametrized as 

poo 

dzif(zi)=i dyif(ci + iyi), (11) 
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Figure 3: The real part of the integrand for the "sunset" diagram in eq. / flOj) . for p 2 = 1, 
m\ = 1 , m| = 4, m 2 = 5. For the left panel, the integration contours have been chosen as 
straight lines parallel to the imaginary axis, and the Z\ integration has already been carried 
out. For the right panel, the contours have been deformed by a rotation in the complex 
plane, corresponding to 9 = 0.4 in eq. i fT3)) . 

where the q are real constants. The integrand contains powers of the kinematic and mass 
invariants as well as gamma functions. 

Let us return to the example of the sunset diagram, eq. (fTUl) . For Euclidean momentum, 
p 2 < 0, the dependence of the invariants on the integration variables yi is given by simple 
oscillating exponentials, e.g. {m\)' iy2 , while the gamma functions rapidly tend to zero when 
their arguments take values away from the real axis. Therefore, in this case, the integrals 
are well behaved for numerical integration, and a moderate finite integration interval G 
[—10, 10] is sufficient to reach a high precision for the result. 

However, for physical momentum, p 2 > 0, an exponentially growing term appears in the 
integrand: 

(-p 2 ) z:i = (p 2 ) C3+iy3 (-l - ie) C3+iy:i = (p 2 )<*+*va e -**<* e *va m ( 12 ) 

Thus the integrand is an oscillating function with slowly decreasing amplitude for y$ — > oo. 
While the integral is formally convergent, numerical integration routines have difficulties 
dealing with this situation, since one needs to integrate to very large values of y%, i. e. over 
many oscillation periods, to reach acceptable precision. See Fig. [3] for illustration. 

The problem of slowly converging oscillating tails of the integrand is generic for any MB 
representation of a loop integral with physical (Minkowski) momenta. Special techniques 
for numerical integration of oscillating integrands like Filon's method [19] fail since the 
oscillation period varies substantially within the integration range. Instead, as shown in 
the following, suitable variable transformations can substantially improve the convergence 
properties. 

Next, one can exploit the freedom to deform the integration contours in the complex 
plane, as long as no contour crosses a pole of the integrand. The latter condition is rather 
restrictive, but if all contours are rotated in the complex plane by the same angle it is 
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guaranteed that no pole crossing occurs: 

q + iyi -> q + (9 + i)yi. (13) 

In eq. f lT2"j) this leads to 

so that the problematic exponential exp(ny 3 ) can, in principle, be compensated by the a 
suitable choice of 9. For the example shown in Fig. [3] (right), it is evident that with this 
contour deformation the integrand drops off very fast away from the origin, and thus is much 
better suited for numerical integration. 

For a mult i- dimensional MB integral, the choice of 9 is more delicate since this one 
parameter affects all integration variables. It helps to replace the variables yi t ,.. tn by hyper- 
spherical coordinates, i.e. n — 1 angles and one radial coordinate: 

Vi = r cos 0i, 

y 2 = r sin 0i cos 02, 



(15) 



y n _i = r sm 0i • • • sin 0„_ 2 cos <p n _ u 
y n = r sin 0i ■ • • sin n _ 2 sin n _i. 



In this representation, only r is affected by the modification (TIB"]) of the integration contours, 
so that one can in principle choose 9 such that all problematic exponentials of the form (1121) 
are cancelled, function of the 0,-. 

After applying the variable transformations, the multi-dimensional MB integral is suffi- 
ciently well-behaved that it can be evaluated with standard Monte-Carlo integration algo- 
rithms like Vegas [2oJ3- In particular, the Monte-Carlo error decreases systematically with 
the number of integration points, which is an important requirement for the convergence of 
the numerical integration routine. 

However, MB integrals with high dimensionality can still require a large amount of com- 
puting time to reach a result with acceptable precision. For two-loop diagrams one can 
obtain MB integrals with more than 10 dimensions, depending on the topology and the 
number of independent masses and momenta involved. Fortunately, it turns out that in 
most cases some of the integrations can be performed analytically. For this purpose, the 
convolution theorem for Mellin transforms proves to be useful. The Mellin transform of a 
function / is defined as 

/■oo 

M[f(x);z]= dxf{x)x z -\ (16) 







For the Mellin transform of a product of two functions, the convolution theorem then states 

i rc+ioo 

M[f(x)g(x);s} = — / dzM[f{x) ] z]M[g{x);s-z], (17) 

Z7TI J c—ioo 

*For relatively simple integrals with four or less integration dimensions, deterministic integration methods 
like Gaussian quadrature prove more efficient. 
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where the integration is parallel to the imaginary axis and the integrand in eq. (IT7|) must be 
analytic for Re(z) in some region around c. 

With suitable choices of / and g one can derive formulas for MB integrals of two or four 
gamma functions: 
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27T? J c—ioo 

= r(a)r % r ( ( ^ + 7 " a) ^(°.7. + 7, 1 - «) [0 > c > 1 - Re(«, 1 - Re( 7) ] 

where in each case the contours are assumed to be straight lines parallel to the imaginary 
axis satisfying the conditions in square brackets, and 2-P1 is the Gauss hypergeometric func- 
tion. These formulas should be applied before the variable transformations. The maximal 
reduction of the integration dimension is achieved if the formulas are used both before and 
after the extraction of the 1/e poles, since in each case one obtains different combinations 
of terms that match the left-hand sides of eqs. (|T8|) - (!20|) . The reduction procedure has been 
implemented in a computer algebra program based on Mathematica, which automatically 
looks for all possible applications of these equations. 

4 Examples 

In this section numerical results for some examples of scalar two-loop integrals are presented. 
The generation of the MB representation, the extraction of the 1/e poles, and the application 
of the reduction formulas, as described in the previous section, has been implemented in 
Mathematica. For the numerical integration a C code has been developed, which uses 
Gaussian quadrature for low dimensionality and the Vegas Monte-Carlo routine for high 
dimensionality. 

In order to test the method and its computer implementation, results will be shown first 
for several known two-loop integrals, which have been calculated earlier with other methods. 
A few examples that cannot be computed with any previously published method will be 
presented at the end of this section. 

It should be stressed that all examples in this section cannot be evaluated numerically 
with MB representations by using the standard representation eq. (TTTj) . as implemented e. g. 
in the MB package of Ref. [T7], see also Ref. [TB]. The reason is that for physical momenta, 
Pi > 0, the integrands have long oscillating tails that converge too slowly. Therefore, for 
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(p 2 ,m?,m|,m§) 


(1,1,4, 5) 


(16, 2, 2, 1) 


Result from MB repr. 
(with integration error) 


18.19558(2) 


17.5083(1) - 0.1075(l)i 


Result from disp. rel. 


18.19559 


17.50812 - 0.10747i 



Table 1: Numerical results for the finite part of the sunset diagram in Fig. \2i obtained 
with the MB representation method described in this report (second line) and the dispersion 
relation method of Ref. (third line). Monte- Carlo integration errors for the last quoted 
digit are given in parentheses. The scale of dimensional regularization has been chosen fi = 1. 



all cases shown in the following, the variable transformations and contour deformations 
introduced in the previous sections have been used. 



Self-energy topologies: The scalar sunset diagram in Fig. [2] leads to a three-dimensional 
MB integral, as shown in section 121 together with one- and two-dimensional integrals origi- 
nating from pole crossings when the limit e — > is taken. Expanding in powers of e, it turns 
out that the UV poles are given by terms without integrals: 



ml + m\ + m\ 1 
2F 2 + e 



m tU-lE - logmj - — 



2,fin" 



(21) 



Numerical results for the finite part are shown in table [T] for two choices of parameters, 
comparing with values obtained through the dispersion relation method of Ref. |12j . 

The scalar five-propagator self-energy in Fig. H] can be cast into a ten-dimensional MB 
representation, which can be reduced to integrals with dimension up to six with the help of 
the reduction formulas presented in the previous section. Note that this integral is UV-finite, 
so that no terms from pole crossings occur. 

Owing to the relative large dimensionality of the integration, the numerical integration 
is much slower than in the previous example. With about one hour evaluation time on a 
single core of a Pentium IV computer with 2.4 GHz a precision of 0.1% to 1% is obtained, 
depending on the values of the masses and the momentum, see table [2J They compare well 
with results based on the dispersion relation method of Ref. [T2] . 



Vertex topologies: Figure [5] shows a few examples for two- loop vertex contributions with 
three or four different scales. The three scalar diagrams lead to MB integrals with 9, 10, 
and 11 dimensions, respectively, which can be reduced to five dimensions each using the 
reduction formulas in eqs. (fl8|) ff. In table [3j numerical results are compared with numbers 
obtained with the dispersion relation method of Refs. 0H2] . 

Box topologies: This subsection will discuss some two-loop box integrals corresponding 
to typical topologies for next-to-next-to-leading order QCD corrections to pp — > ti, shown in 
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Figure 4: Five- propagator two-loop diagram. As before, all propagators are assumed to be 
scalars. 



(p 2 , ml, mj, ml, ml, ml) 


(-1, 1, 1, 1, 1, 1) 


(1, 3.7, 3.7, 1, 3.7, 3.7) 


(20, 1, 1, 1, 1, 1) 


Result from MB repr. 
(with integration error) 


-0.6809(4) 


-0.246(2) 


0.439(1) + 0.259(7)i 


Result from disp. rel. 


-0.68088 


-0.25421 


0.43952 + 0.25365z 



Table 2: Numerical results for the self-energy diagram in Fig. dj obtained with the MB 
representation method described in this report (second line) and the dispersion relation 
method of Ref. jJ2| (third line). Monte-Carlo integration errors for the last quoted digit are 
given in parentheses. The scale of dimensional regularization has been chosen fi = 1. 




(a) (b) (c) 

Figure 5: Planar two-loop vertex diagrams which appear as part of the two-loop corrections 
to Z —>■ bb. As before, all propagators are assumed to be scalars, and external legs are 
constrained to be on-shell. 



Diagram 


Fig. E (a) 


Fig. E (b) 


Fig. E (c) 


Result from MB repr. 
(with integration error) 


(2.24(2) + 2.34(2)i) x 10" 9 


1.04(1) x 10~ 9 


3.78(3) x 10~ 10 


Result from disp. rel. 


(2.241 + 2.376i) x 10" 9 


1.036 x 10" 9 


3.805 x 10~ 10 



Table 3: Numerical results for the vertex diagrams in Fig. El obtained with the MB repre- 
sentation method described in this report (second line) and the dispersion relation method 
of Ref. (third line). Monte-Carlo integration errors for the last quoted digit are given 
in parentheses. The following values for the masses and regularization scale have been used: 
mz=91, mw=80, m t =175, mH=100, fi=l. 
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(a) (b) (c) (d) 

Figure 6: Some scalar master integrals relevant for two-loop QCD corrections to pp — > 
ti. Thick lines indicate massive top-quark propagators, while thin lines are massless, and 
external legs are constrained to be on-shell. 



(s,t,m 2 ) 


Result from MB representation xlO 15 (with integration error) 


(-100 2 , -100 2 , 175 2 ) 
(400 2 , -100 2 , 175 2 ) 


-2.95397£~ 2 + [50.15(1) + 0.00(l)i]£" 1 - [515.1(4) - 0.0(4)i] 

[17.9514 - 13.6206i]e~ 2 + [4720.1(5) - 668.5(4)i]e _1 
+ [3651(13) - 1481(13)z] 



Table 4: Numerical results for the box diagram in Fig. [6| (c), obtained with the MB 
representation method. Monte-Carlo integration errors are given in parentheses. These 
numbers are based on m t = 1 75 for the top mass and u. = 1 for the regularization scale. 



Fig. [6l In the diagrams massless (thin) lines correspond to light quark or gluon propagators, 
while massive (thick) lines indicate top quark propagators. 

The master integrals corresponding to (a) and (b) in the figure lead to MB integrals with 
dimension up to two and three, respectively, after expansion in e. Therefore they can be 
evaluated very efficiently with deterministic integration methods, like Gaussian quadrature. 
These master integrals were computed previously in Ref. [10J with the method of differential 
equations. Comparing results between the two methodaj, agreement to more than nine 
significant digits has been obtained!- This comparison was carried out both for an unphysical 
parameter point with Euclidean momentum s < 0, as well as for a physical parameter point 
with s > (2m t ) 2 . For these cases, the evaluation of the MB integrals takes only several 
seconds to minutes on a single-core Pentium IV computer. 

To the best of our knowledge, the two-loop integral corresponding to Fig. [6] (c) and (d) 
have not been computed before. Fig.[6](c) leads to a MB representation with a 13- dimensional 
integral, which can be reduced to six dimensions with the help of the reduction formulas. 
Numerical results for a few parameter points are shown in table HI Since the structure of the 
integral is rather complicated, the evaluation is much slower than the previous examples. On 
a single-core Pentium IV computer with 2.4 GHz, the numerical evaluation of one parameter 
point takes 4-5 hours to reach the precision reported in the table. 

Fig. M (d) corresponds to a 9-dimensional MB integral, which decomposes into integrals of 
dimension up to five after application of the reduction formulas and expansion in e. Table [5] 
shows numerical results which have been obtained in about 10 hours evaluation time. 

1"We are grateful to the authors of Ref. [10] for kindly providing numerical results of their calculation. 
* Higher precision could be achieved with more computing time. 
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(s, t, m 2 ) 


Result from MB representation xlO 15 (with integration error) 


(-100 2 , -100 2 , 175 2 ) 
(400 2 , -100 2 , 175 2 ) 


75.8804e- 3 - [1141.55(5) + 0.00(5)i]£- 2 

+ [13338.5(1.4) - 0.5(1.4)i] £ - 1 - [81900(530) + 40(530)i] 

[4.17615 + 12.4156i]e- 3 - [58.04(3) + 272.30(3)i]£- 2 

+ [252.3(9) + 2932.3(9)i]£~ 1 + [742(44) - 20649(44)i] 



Table 5: Numerical results for the box diagram in Fig. \B(d), for the same conditions as in 
Tab.d 



5 Summary 

Mellin-Barnes representations are a powerful tool for the calculation of multi-loop integrals. 
They have been shown to be very useful for systematically extracting poles of the form 
1/(4 — D) in dimensional regularization. The coefficients of these poles and the finite part can 
then be evaluated as numerical integrals. However, for diagrams with massive propagators 
and physical branch cuts, the integrands a priori can be highly oscillatory, so that standard 
numerical integration routines do not yield convergent results. 

In this article, it is shown how the numerical convergence behavior can be improved sub- 
stantially by transforming the integration variables, deforming the integration contours in the 
complex plane, and using analytical formulas to reduce the dimensionality of the integration. 
The usefulness of those modifications was demonstrated by calculating several examples of 
two-loop self-energy, vertex and box topologies, some of which have been computed for the 
first time. 

Several previously known techniques permit much faster evaluation of relatively simple 
multi-loop integrals with few independent masses and momenta than the method presented 
here. However, Mellin-Barnes representations can be applied in principle to any loop integral, 
irrespective of its complexity. With sufficient computing time, arbitrarily precise results can 
be obtained. 

In this work, tensor integrals with external uncontracted Lorentz indices were not dis- 
cussed separately. In fact, for most applications, it is convenient to decompose the loop 
amplitudes into standard matrix elements or form factors, the coefficients of which only con- 
tain scalar integrals. The integrands of these scalar integrals can contain scalar products in 
the numerator involving on one or more of the loop momenta. Such integrals can be evalu- 
ated straightforwardly with Mellin-Barnes representations by expressing the scalar products 
as propagators with negative powers. 
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